##########################################################################################

library(dplyr)
library(data.table)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--maf_path"), type = "character") ,
    make_option(c("--class_type"), type = "character") ,
    make_option(c("--images_path"), type = "character") ,
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--info_public_file"), type = "character") ,
    make_option(c("--class_order_file"), type = "character") ,
    make_option(c("--im_list"), type = "character"),
    make_option(c("--igc_list"), type = "character"),
    make_option(c("--dgc_list"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic"
    maf_path <- paste(work_dir,"/","maf_public",sep="")
	images_path <- paste(work_dir,"/","images",sep="")
	info_file <- paste(work_dir,"/config/tumor_normal.class.list",sep="")
	info_public_file <- paste(work_dir,"/public_ref/combine/MutationInfo.combine.addMolecularSubType.tsv",sep="")
	class_order_file <- paste(work_dir,"/config/Class_order.list",sep="")
	#smg_list <- paste(work_dir,"/public_ref/SMG_sort.list",sep="")

	im_list <- paste(work_dir,"/mutsig_check/im_smg.list",sep="")
	igc_list <- paste(work_dir,"/mutsig_check/igc_smg.list",sep="")
	dgc_list <- paste(work_dir,"/mutsig_check/dgc_smg.list",sep="")
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

maf_path <- opt$maf_path
info_file <- opt$info_file
info_public_file <- opt$info_public_file
images_path <- opt$images_path
class_order_file <- opt$class_order_file
im_list <- opt$im_list
igc_list <- opt$igc_list
class_type <- opt$class_type
dgc_list <- opt$dgc_list

dir.create( images_path , recursive = T )

###########################################################################################

im_maf_file <- paste0(maf_path , "/All_use.IM.maf")
gc_maf_file <- paste0(maf_path , "/All_use.maf")

###########################################################################################
#smg <- data.frame(fread(smg_list , header = F))
dat_im_smg <- data.frame(fread(im_list , header = T))
dat_igc_smg <- data.frame(fread(igc_list , header = T))
dat_dgc_smg <- data.frame(fread(dgc_list , header = T))

class_order <- data.frame(fread(class_order_file , header = T))
info <- data.frame(fread(info_public_file))

dat_im <- data.frame(fread(im_maf_file))
dat_gc <- data.frame(fread(gc_maf_file))

###########################################################################################

smg <- data.frame(Gene_Symbol = unique(c(dat_im_smg$Gene_Symbol , dat_igc_smg$Gene_Symbol , dat_dgc_smg$Gene_Symbol)))

###########################################################################################
## 去除MSI样本
info <- subset( info , !(Molecular.subtype %in% c("POLE" , "MSI") ))
info$Molecular.subtype <- ifelse( info$Molecular.subtype %in% c("EBV" , "unknown") , "Other" , info$Molecular.subtype )
## 去除IGC和DGC都有的样本
info <- subset( info , Class != "IGC + DGC" )

###########################################################################################
## 一个人多个样本算一个样本
dat_im$Tumor_Sample_Barcode <- paste0( dat_im$Tumor_Sample_Barcode , "_IM" )

info_im <- subset( info , From == "NJMU"  )
info_im$Tumor <- paste0( info_im$Tumor , "_IM" )
info_im$Class <- "IM"

info <- rbind(info_im , info)

info$Class <- factor(info$Class , levels = c("IM", "IGC" , "DGC" ) , ordered=T)
info$Normal <- gsub( "_IM" , "" , info$Tumor )
rownames(info) <- info$Tumor

###########################################################################################

CreateMutMatrix <- function(dat = dat , smg = smg , Variant_Type = Variant_Type ){

	mut <- dat
	mut$Variant_Classification <- as.character(mut$Variant_Classification)
	mut <- mut[which(mut$Variant_Classification %in% Variant_Type),]

	## INDEL 和 INS 合并
	mut[grep("In_Frame",mut$Variant_Classification),'Variant_Classification'] = "In_Frame"
	mut[grep("Frame_Shift",mut$Variant_Classification),'Variant_Classification'] = "Frame_Shift"

	## Gene只考虑有突变的SMG和CGC
	Gene <- smg[smg$Gene_Symbol %in% mut$Hugo_Symbol,]

	## 构建矩阵
	Sample <- info[,"Tumor"]
	maf_matrix <- matrix("" , ncol = length(Sample) , nrow = length(Gene) , dimnames = list(Gene,Sample))

	for(gene in rownames(maf_matrix)){
		print(gene)

		## 突变类型
		for(tumor in colnames(maf_matrix)){
			index <- which(mut$Hugo_Symbol==gene & mut$Tumor==tumor)
			if(length(index)==0){
				var=""
			}else if(length(index)==1){
				var <- mut[index,'Variant_Classification']
			}else if(length(index)>1){
				var <- "Multiple_Hits"
			}
			print(var)
			maf_matrix[ which(rownames(maf_matrix)==gene) , which(colnames(maf_matrix)==tumor) ] <-  var
		}
	}

	return(maf_matrix)
}

###########################################################################################
## 产生输入矩阵
Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
Variant_Type_Combine <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift","In_Frame","Splice_Site","Nonstop_Mutation","Multiple_Hits")

###########################################################################################
## 合并所有的MAF文件
dat <- rbind( dat_gc , dat_im )
dat <- dat[,colnames(dat) != "Variant_Type" ]
dat$Tumor <- dat$Tumor_Sample_Barcode

MutMatrix <- CreateMutMatrix( dat = dat , smg = smg , Variant_Type = Variant_Type )

###########################################################################################

use_sample <- subset( info , Class==class_type )$Tumor
MutMatrix <- MutMatrix[,use_sample]

##############################################################################
result_dat <- c()
for( i in 1:(nrow(MutMatrix)-1) ){

	gene_1 <- rownames(MutMatrix)[i]

	for( j in (i+1):nrow(MutMatrix) ){
		gene_2 <- rownames(MutMatrix)[j]

		## 两者均有
		a <- length(which(MutMatrix[i,]!="" & MutMatrix[j,]!=""))
		## 只有gene_1有
		b <- length(which(MutMatrix[i,]!="" & MutMatrix[j,]==""))
		## 只有gene_2有
		c <- length(which(MutMatrix[i,]=="" & MutMatrix[j,]!=""))
		## 都没有
		d <- length(which(MutMatrix[i,]=="" & MutMatrix[j,]==""))

		result=fisher.test(matrix(c(a,b,c,d),nrow=2))

    p=result[["p.value"]]
    OR=round(result[["estimate"]][["odds ratio"]],3)

    tmp <- data.frame( geneA = gene_1 , geneB = gene_2 , ConMut = a , Only_geneA_Mut = b , Only_geneB_Mut = c , NoMut = d , p = p , OR = OR )

    result_dat <- rbind( result_dat , tmp )
	}
}

images_name <- paste0(images_path , "/MutuallyExclusive.",class_type,".tsv")
write.table( result_dat , images_name , row.names = F , quote = F , sep = "\t" )